#include <QtCore/QCoreApplication>
#include <QPluginLoader>
#include <QDir>
#include <QFile>

#include "../plugin.h"
#include "../data/fileplugin.h"
#include "../data/reg3dvolumedata.h"
#include "../contree/contreeconstructor3d.h"
#include "../contree/augmentedctree.h"
#include "../contree/ctreecontourextractor.h"
#include "../math/timer.h"
#include "../mact/dctconstructor.h"
#include "../mact/dualgraph.h"
#include "../shape/sphere.h"
#include "../contour/mcconextractor.h"
#include "../sdf/geometryflow.h"
#include "../sdf/distancetransform.h"
#include "../sdf/pocketextractor.h"
#include "../contour/molskin.h"

#include <iostream>
#include <vector>
#include <algorithm>
#include <utility>

using namespace PBVP;

class MaximumPoint {
public:
	MaximumPoint(float _x = .0f, float _y = .0f, float _z = .0f, float _f = .0f) {
		x = _x;
		y = _y;
		z = _z;
		f = _f;
	}

	float x, y, z;
	float f;

	friend bool operator < (const MaximumPoint& p1, const MaximumPoint& p2) {
		return p1.f < p2.f;
	}
};


int main(int argc, char *argv[])
{

	QCoreApplication a(argc, argv);
	//return a.exec();

	// load plugins
	QDir pluginsDir = QDir(a.applicationDirPath());
	std::cout << pluginsDir.absolutePath().toStdString() << std::endl;
	if (pluginsDir.dirName().toLower() == "debug"
		|| pluginsDir.dirName().toLower() == "release") {
			pluginsDir.cdUp();
	}
	// pluginsDir.cdUp();
	pluginsDir.cd("plugins");

	QStringList pluginEntries = pluginsDir.entryList(QDir::Files);

	FilePluginInterface* rawivPlugin;
	FilePluginInterface* rawPlugin;
	FilePluginInterface* rawncPlugin;
	FilePluginInterface* pdbPlugin;

	// check additional potential plugin path if necessary
	if(pluginEntries.size() == 0) {
		pluginsDir = QDir(a.applicationDirPath());
		pluginsDir.cd("plugins");
		pluginEntries = pluginsDir.entryList(QDir::Files);
	}

	std::cout << "Loading plugins in dir" << pluginsDir.absolutePath().toStdString() << std::endl;

	foreach (QString fileName, pluginEntries) {
		if (fileName.contains("rawiv"))
		{
			QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
			QObject *plugin = loader.instance();
			if (plugin) {
				rawivPlugin = qobject_cast<FilePluginInterface *>(plugin);
				if (rawivPlugin) {
					rawivPlugin->initialize(NULL);
					std::cout << rawivPlugin->description().toStdString() << std::endl;
				}
			}
		} else if (fileName.contains("rawfile"))
		{
			QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
			QObject *plugin = loader.instance();
			if (plugin) {
				rawPlugin = qobject_cast<FilePluginInterface *>(plugin);
				if (rawPlugin) {
					rawPlugin->initialize(NULL);
					std::cout << rawPlugin->description().toStdString() << std::endl;
				}
			}
		} else if (fileName.contains("rawncfile"))
		{
			QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
			QObject *plugin = loader.instance();
			if (plugin)
			{
				rawncPlugin = qobject_cast<FilePluginInterface *>(plugin);
				if (rawncPlugin)
				{
					rawncPlugin->initialize(NULL);
					std::cout << rawncPlugin->description().toStdString() << std::endl;
				}
			}
		} else if (fileName.contains("pdbfile"))
		{
			QPluginLoader loader(pluginsDir.absoluteFilePath(fileName));
			QObject *plugin = loader.instance();
			if (plugin)
			{
				pdbPlugin = qobject_cast<FilePluginInterface *>(plugin);
				if (pdbPlugin)
				{
					pdbPlugin->initialize(NULL);
					std::cout << pdbPlugin->description().toStdString() << std::endl;
				}
			}
		}
	}

	char fname[256];
	Timer timer;
	timer.getCurrentTime();
	std::cout << "processing file " << argv[1] << std::endl;
	ProteinData* pdb = (ProteinData *) (pdbPlugin->loadFile(argv[1]));
	QVector<Atom> atoms = pdb->getAtoms();
	MolSkin molskin;
	GeometryData* pMolSurf = molskin.constructMolSurface(atoms);
	sprintf(fname, "%s_ms.raw\0", argv[1]);
	rawPlugin->saveFile(pMolSurf, fname);

	PocketExtractor pocketExtractor(40, 1.5, 0);		// Don't filter
	Reg3DVolumeData *pPocVol = pocketExtractor.computePocketVolume(pMolSurf);
	sprintf(fname, "%s_pks.rawiv\0", argv[1]);
	rawivPlugin->saveFile(pPocVol, fname);
	delete pMolSurf;
	delete pdb;

	AugmentedContourTree *acTree = new AugmentedContourTree();

	sprintf(fname, "%s_pks.rawiv.ctr\0", argv[1]);
	if (!acTree->loadFile(fname))
	{
		ConTreeConstructorReg3D *construcCT = new ConTreeConstructorReg3D(pPocVol);
		delete acTree;
		acTree = construcCT->getAugmentedContourTree(); 
		delete construcCT;
		acTree->saveFile(fname);
	}
	CTreeContourExtractor* pConExtractor = new CTreeContourExtractor(pPocVol, acTree);
	pConExtractor->setInnerRegion(1);
	std::vector<GeometryData*> contours;
	int nCon = pConExtractor->extractLargestContours(2.0, contours, 4);		// try to extract the four largest pockets
	delete pConExtractor;
	delete acTree;

	unsigned int dim[3];
	pPocVol->getDimension(dim);
	
	for (int i = 0; i < dim[0]*dim[1]*dim[2]; i++)
	{
		float val = pPocVol->getValue(i)*pPocVol->getValue(i);
		// float val = fabs(pReg3Data->getValue(i));
		pPocVol->setValue(val, i);
	}
	
	
	MCContourExtractor conExtractor;
	int nflow = 100;
	for (int i = 0; i < nCon; i++)
	{
		sprintf(fname, "%s_pk%d.raw\0", argv[1], i+1);
		rawPlugin->saveFile(contours[i], fname);
		
		DistanceTransform dst(contours[i]);
		dst.setDimension(dim);
		double minExt[3], maxExt[3];
		pPocVol->getMinMaxExt(minExt, maxExt);
		dst.setMinMaxExt(minExt, maxExt);
		Reg3DVolumeData* pShape = dst.computeSDF();

		for (int v = 0; v < dim[0]*dim[1]*dim[2]; v++)
		{
			pShape->setValue(-pShape->getValue(v), v);
		}
		GeometryFlow flow(pShape, pPocVol);	
		flow.evolve(nflow);
		GeometryData *pNewPocket = conExtractor.extractContour(pShape, 0);
		sprintf(fname, "%s_pk%d_evolved%d.raw\0", argv[1], i+1, nflow);
		rawPlugin->saveFile(pNewPocket, fname);
		delete pNewPocket;

		delete pShape;
	}

/*
	Reg3DVolumeData *pReg3Data = (Reg3DVolumeData *)(rawivPlugin->loadFile(argv[1]));

	// look for maxima in the volume function
	// For pocket function, we want f > 0 and no other maximum in
	// a sphere of radius f
	std::vector<MaximumPoint> maxvec;
	std::vector<MaximumPoint>::iterator it;
	unsigned int dim[3];
	pReg3Data->getDimension(dim);
	for(int k = 0; k < dim[2]; k++) {
		for(int j = 0; j < dim[1]; j++) {
			for(int i = 0; i < dim[0]; i++) {
				float f = pReg3Data->getValue(i, j, k);
				bool maxflag = true;
				for(int kk = std::max<int>(0, k-1); kk < std::min<int>(k+1, dim[2]); kk++) {
					for(int jj = std::max<int>(0, j-1); jj < std::min<int>(j+1, dim[1]); jj++) {
						for(int ii = std::max<int>(0, i-1); ii < std::min<int>(i+1, dim[0]); ii++) {
							if(ii != i || jj != j || kk != k) {
								float f2 = pReg3Data->getValue(ii, jj, kk);
								if(f2 >= f) maxflag = false;
							}	// if statement
						}
					}
				}
				if(maxflag && f > 0.5) {	// consider maximum points with SDF > 0.5
					float coord[3];
					pReg3Data->getCoordinates(coord, i, j, k);
					MaximumPoint mp(coord[0], coord[1], coord[2], f);
					maxvec.push_back(mp);

				}
			} // for loops
		}
	}
	std::sort(maxvec.begin(), maxvec.end());
	//for(int i = 0; i < maxvec.size(); i++) {
	//	std::cout << "MP #" << i << ": " << maxvec[i].x << " " << maxvec[i].y << " " << maxvec[i].z << " " << maxvec[i].f << std::endl;
	//}
	int count = 0;
	GeometryData *pSpheres = new GeometryData();
	int size = maxvec.size();


	// consider the pocket containing the largest sphere now
	MaximumPoint mp = maxvec.back();	// the largest sphere
	Vector3 v(mp.x, mp.y, mp.z);
	std::vector<MaximumPoint> pocket1;
	pocket1.push_back(mp);
	for (int i = 0; i < size; i++)
	{
		MaximumPoint p2 = maxvec[i];
		Vector3 v2(p2.x, p2.y, p2.z);
		if(v.distance(v2) < (mp.f + p2.f)) {
			pocket1.push_back(p2);
			//Sphere *s = Sphere::makeSphere(v2, p2.f);
			//GeometryData *pGeo = s->toGeometry(0);
			//pSpheres->merge(pGeo);
			//delete pGeo;
			//delete s;
		}
	}
	// construct the union of the spheres in the largest pocket
	float orig[3], span[3];
	Reg3DVolumeData *pSphVol;
	pReg3Data->getOrigin(orig);
	pReg3Data->getSpan(span);
	printf("span = (%f %f %f)\n", span[0], span[1], span[2]);


	GeometryData *pPocket =  (GeometryData *)(rawPlugin->loadFile("1HOS_pk1.raw"));
	DistanceTransform dst(pPocket);
	dst.setDimension(dim);
	double minExt[3], maxExt[3];
	pReg3Data->getMinMaxExt(minExt, maxExt);
	dst.setMinMaxExt(minExt, maxExt);
	pSphVol = dst.computeSDF();
	for (int i = 0; i < dim[0]*dim[1]*dim[2]; i++)
	{
		pSphVol->setValue(-pSphVol->getValue(i), i);
	}
	// compute the zero contour of the volume function of spheres
	MCContourExtractor conExtractor;
	pSpheres = conExtractor.extractContour(pSphVol, 0.2);
	rawPlugin->saveFile(pSpheres, "spheres.raw");
	rawivPlugin->saveFile(pSphVol, "spheres.rawiv");
	delete pSpheres;
	// square the pocket SDF
	int nxyz = dim[0]*dim[1]*dim[2];
	for (int i = 0; i < nxyz; i++)
	{
		float val = pReg3Data->getValue(i)*pReg3Data->getValue(i);
		// float val = fabs(pReg3Data->getValue(i));
		pReg3Data->setValue(val, i);
	}
	GeometryFlow flow(pSphVol, pReg3Data);
	for (int i = 0; i < 50; i++)
	{
		flow.evolve(1);
		printf("Iteration # %d\n", i+1);

		pSpheres = conExtractor.extractContour(pSphVol, 0);

		sprintf(fname, "spheres%d.raw\0", i+1);
		rawPlugin->saveFile(pSpheres, fname);
		// sprintf(fname, "spheres%d.rawiv\0", i+1);
		// rawivPlugin->saveFile(pSphVol, fname);

		if (i > 0 && i%10 == 0)
		{
			flow.reinitilization();
		}
		delete pSpheres;
	}


	std::cout << "Elapsed time = " << timer.getElapsedTime() << std::endl;
	delete pReg3Data; */
}
